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Abstract 

This work concentrates on numerical simulations of Photonic Crystal structures using 
basis-expansion methods, with a main focus on simulating disorder. The plane-wave and 
guided-mode expansions are outlined as tools to compute the Bloch modes of a structure, on 
the basis of which the Bloch-mode expansion formalism is outlined - the latter allowing for 
simulations of large structures in presence of disorder. As a first illustration of the method, we 
apply it to three gentle-confinement cavities, to obtain results for their quality factors similar 
to the theoretically predicted in the literature. Furthermore, we compute that random disorder 
can drive those factors down to the experimentally measured values. As a second application, 
we study the effect of irregular hole shapes in a PHC waveguide, and find that the correlation 
length of the irregularity (i.e. the typical scale of the roughness of the features) matters: 
for higher correlation lengths, the computed modes show both higher band broadening and 
higher loss rates. 
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Introduction 


1.1 What is a photonic crystal? 

It has been shortly less than a quarter of a century since what are now considered some of the 
seminal works on photonic crystals, n, m first appeared, but the signihcance of those structures 
was so quickly realized that the field is currently in the spotlight of photonics research and de¬ 
velopment. A photonic crystal, in essence, is any form of a medium with periodically varying 
dielectric constant, but more generally it can be a structure with such an underlying periodicity, 
which is however broken by e.g. unwanted disorder or deliberately placed defects for the achieve¬ 
ment of particular tasks. Such crystals also appear in nature, but it is the ability to design and 
manufacture them in the lab which shows a potential for numerous applications in optical devices. 

A general classification that can be made between different photonic crystals is the dimen¬ 
sionality in which the periodicity occurs; namely, there can be ID photonic crystals, where the 
periodicity is in one dimension only, and the crystal is homogeneous in the other two dimensions - 
an example for this is a collection of alternating layers of two different materials. Similarly, a 2D 
crystal would be periodic in two spacial directions and homogeneous in the third - an example can 
be periodically occurring dielectric rods in air (fig. (a)). And finally, a 3D crystal is described 
by some periodic dielectric constant in all directions - for example tangential dielectric spheres in 
air (hg. [l](b)). 



Figure 1: (a): GaAs/AlGaAs pillars in air, image from [3]. (b): 3D photonic crystal with the 
“opal structure” of dielectric spheres in air, image from [3]. (c): waveguide created by a “missing” 
row of air holes in a dielectric slab, image from [5]. 


Photonic crystals (PHCs) have many similarities to the standard crystals studied in Solid State 
Physics. Generally, the periodicity of the dielectric function can be represented in some form as a 
lattice of primitive cells. Photons in a PHC thus behave in a similar way to electrons in a crystal 
lattice, both particles being governed by an underlying wave equation. We can thus construct 
band diagrams for a PHC, just as we do for a standard crystal. The periodicity of the structure 
of the PHC often provides a band-gap over a given frequency range which provides a basis for 
interesting properties and applications. The main distinction between the two types of crystals 
is that, first, electrons are fermions, while photons are bosons, so they essentially obey different 
statistics, and second, that photons do not carry electric charge, so we do not need to consider 
Coulomb interaction, as opposed to when considering electrons and holes. 

1.2 Current developments and applications 

For a detailed outline of photonic crystal development and application possibilities, refer to Chap¬ 
ter 1.4 of [5], and references therein; here we will briefly highlight the areas of science where PHCs 
can prove useful. 
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• Spontaneous radiation management: the main topic of [2] and more recently of a detailed 
review paper by S. Noda [7], this is arguably the area where PHCs will prove most useful. 
In a 3D PHC, there can be a band gap in all spatial directions such that control as strong as 
complete inhibition of spontaneous emission of a material introduced into such a structure 
can be achieved. However, even 2D structures present interesting opportunities for control 
of spontaneous emission, e.g. when the emission is inhibited in two of the spatial dimensions 
so that all radiation energy is redirected to only exit through the third. An overview of 
the applications of the use of PHCs for spontaneous emission control includes improving 
semiconductor lasers or solar cell efficiency, and optical cavities based on PHCs can be a 
basis for new types of radiation sources for illumination, display and optical communication 
purposes. 

• Introducing defects in a PHC provides numerous possibilities for constructing waveguides 
(fig- a (c)), microcavities, splitters, couplers and combiners. The waveguides are usually 
produced by introducing a line defect, and they work due to the presence of the photonic 
band gap, rather than due to internal reflection. This allows them to form very sharp bends, 
e.g. at an angle 7r/2. A splitter is essentially a number of such waveguides connected at 
a single point. Since we have great freedom in varying the waveguide parameters, we can 
finely control the operation of the splitter. 

• PHCs can also be used as super-prisms, super-lenses, multiplexers and demultiplexers for 
dispersion management. 

• Nonlinear optics applications: for example for optical storage elements or logical elements. 

• Slow light applications: the group velocity at a specific wavelength in a PHC can be very 
low. 

• Optical fibers: PHCs provide many possibilities to improve existing technology, or to provide 
a basis for a completely new type of fibers, where light is localized due to the photonic band 
gap and an introduced defect. 

• Quantum computing: PHCs can be potentially valuable for manufacturing quantum compu¬ 
tation elements e.g. for linear optics quantum computation as proposed in [8]. In addition, 
at least theoretically, quantum dots used as qubits can be embedded in a PHC structure, 
and coupling between individual dots (e.g. in order to perform two-qubit operations) can 
be achieved through the electromagnetic modes of the underlying structure, which can be 
very finely tuned. As is the case in virtually every area of development related to quantum 
computers, there is a myriad of technical challenges in realizing such a setup in practice, but 
it is definitely worth the effort as it has certain advantages over some of the other qubit and 
qubit gates realizations. 

The core of this thesis is contained in sections 2 and 3: in section 2, we outline some of the 
possibilities to numerically compute the properties of certain PHC structures. Section 3 introduces 
the Bloch-mode expansion method, a tool which is well-known and established in the field of Solid 
State Physics, but has only recently been applied in the PHC case [3]. We then use this method 
to obtain some interesting results regarding disorder (both designed and random) in the crystals. 
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2 Simulating a regular structure 

The starting point for a theoretical treatment of a Photonic Crystal, when no quantum mechanics 
effects need be taken into account, are the Maxwell equations of electrodynamics. As long as we 
are interested only in the free light propagation in the crystal, we can assume zero external charges 
and electrical currents, which allows us to write the Maxwell equations for the electric field E and 
the magnetic field H as 


V • (£E) = 0, 


V H = 

V X E = 

V X H = 


0 , 

c dt ’ 
e^E 
c dt ’ 


( 1 ) 


where e = e(r) is the dielectric permittivity and /i = /i(r) is the magnetic permeability of the 
medium, which in general depend on r. In fact, it is their periodic dependence with r which is in 
the very definition of a PHC structure. Everywhere from now on, we will assume /i(r) = 1, which 
is true for a wide range of materials. If we also restrict the investigation to monochromatic light 
of frequency w, the time dependence will simply be given by a phase shift of wt, and we are left 
with the following equations for the electric and magnetic fields: 


X V X E(r) = ‘^E(r), (2) 

e(r) 

1 

V X X H(r) = ^H(r). 

£(r) 

In general, there are several popular ways to approach numerically the problem of hnding 
the properties of a certain PHC structure. The Finite-Difference Time-Domain method (FDTD) 
is a brute-force solution of Maxwell’s equations over a discretized space- and time-grid; while 
it is a very general method based on a first-principle computation, it is unfortunately often too 
demanding in computational resources. This is partly due to the fact that the discretization 
of space and time can introduce unwanted artifacts, and it is only in the limit of infinitesimal 
intervals that the solution is expected to correspond to the full Maxwell equations solution. In 
the “scattering-matrix” method, the modulation of the dielectric constant is represented as a 
system of scattering sources, repeated periodically in agreement with the PHC periodic structure. 
The scattering-matrix formalism (encountered frequently in quantum mechanics) is then used to 
investigate the propagation of light through the medium. This method, while being a bit more 
involved mathematically than the FDTD, is still essentially a first-principle simulation and so is 
expected to produce very good results, but it can also easily become computationally demanding 
to apply. 

The class of methods this thesis specializes in is the “mode expansion” methods, in which 
the light mode of an unknown system is expanded over the basis of known functions. Such a 
method is not very ubiquitous in the sense that one should choose a proper basis (which is not 
always easy to find) for every structure at hand, and that it involves an infinite summation which 
is not guaranteed to converge for any finite number of terms taken into account. Still, in cases 
where such a method works (and there are some in which it works perfectly), it proves to be 
very advantageous to use it, because it would usually have a lower computational complexity, and 
furthermore it could provide insight into the physics of the system by showing the degree of mixing 
between the actual mode of the system and the basis modes. 

In the rest of this thesis, we will concentrate on a PHC formed by cutting circular holes of 
radius i?, on a lattice (triangular or square) with periodicity a in a dielectric slab of thickness 


3 



d with dielectric constant £ 2 - In general, the outside medium can be thought of as having some 
dielectric constant £ 1 , but in all the actual computations we stick to air and so set £1 = 1. Such 
structures are currently easy to produce in the lab with standard lithographic techniques, and 
thus provide a good research focus. 

2.1 Plane-wave expansion 

In standard solid state structures, the most common basis over which to expand the electron modes 
is the set of plane-wave functions of momentum k, '0k(r) = A straightforward analogy that 
can help us obtain the band diagram of a regular PHC is by plane-wave expansion (PWE) of the 
electromagnetic modes. To do that, we first invoke the Bloch theorem, which can be proven to 
apply to the classical electromagnetic wave-equations in just the same way it applies for particles 
governed by the Schrodinger wave-equation. Its main result is that the modes of the electric and 
the magnetic fields of a periodic structure can be written as 


E(r) = Ek„(r) • (3) 

H(r) = Hk„(r) • 

where the functions Ekn(r) and Hkra(r) are called the Bloch modes of the structure. They 
have periodicity of the underlying lattice, and are labeled by their momentum k and a band index 
n. Based on the periodicity of the lattice, one can also define inverse lattice vectors G such that 
the Bloch modes corresponding to momenta k and k G are equivalent; thus, generally, only the 
modes within an area of the fc-space where they can be uniquely defined (e.g. the first Brillouin 
zone) are considered. Finally, those modes can be represented in Fourier space by a summation 
over all reciprocal lattice vectors G, or in other words they can be expanded on the basis of plane 
waves with momenta k -|- G: 


Ek„(r) =^Ek„(G)e*('^+°)•^ 

G 

Hk„(r)=^Hk„(G)e*('^+°)-'. 

G 


(4) 


2.1.1 Numerical procedure 

Here and in the rest of this thesis, we will consider only 2D-crystals. If in the third dimension 
translation symmetry is kept, i.e. if s{p,z) = e{p,0) Vz and so d = 00, and we restrict the 
discussion to light with in-plane wave-vector, there are two polarization states possible, namely 
Transverse Electric (TE), in which the z-component of the electric field vanishes as do the in-plane 
components of the magnetic field, and Transverse Magnetic (TM), in which the magnetic field is 
in-plane while the electric field is in the z-direction only. For TE-polarization, one obtains from eq. 
(2) the so-called Helmholtz equation for the perpendicular magnetic field component (we assume 
the periodicity of the crystal is in the x-y-direction) [TO] : 


did 

dx £(r) dx 


did 

dy £(r) dy 


Hz{r) = 


(5) 


This equation holds for any mode supported by the system, so in particular - for the Bloch 
modes Hkn- It can be handled more conveniently in Fourier space, where eq. ([^ can be rewritten 
as the eigenvalue problem (often called the “master equation”). 


^£-i(G-G')(k-bG)(k-bG')i?.,kn(G') = (‘^)'i?.,kn(G), (6) 
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where i?z,kn(G) simply denotes the z-component of Hkn(G), which for TE-polarization is the 
only non-vanishing one, and we have used the Fourier expansion of the inverse of the permittivity 
function, e“^(G). Mathematically, this is simply 

(G) = ^ ^ ^ exp(-iGr)d2r, (7) 

with A being the area of the unit cell. In the case of circular holes in a dielectric medium 
of permittivity 62 , the integral in eq. 0 can be computed analytically and represented using a 
Bessel function of the first kind, namely 


e-i(G^0)=2/ 



1 \ Ji(Gr) 
62 ) Gr 


£1 £2 


( 8 ) 


where ei is the permittivity of the surrounding medium (usually air), and / is the so-called 
’’filling factor”, defined as the ratio of the hole area to the total elementary cell area, so in this case, 
/ = !A. However, due to the discontinuities of the permittivity at the edges of the holes, it 

turns out that a method to compute £“^(G —G') which was first used by Ho [TT] gives much better 
numerical convergence. Namely, the method involves computing the Fourier transform e(G — G') 
as a matrix over all possible combinations of differences of G-vectors taken into account. 


/£(Gi-Gi) 

£(G2-Gi) 

£(Gi-G2) • 

£(G2-G2) • 

• £(Gi-G„)\ 

• £{G 2 - G„) 

)^e(G„ — Gi) 

£(G„-G2) • 

• £(G„-G„)/ 


(9) 


and then using matrix inversion to obtain £“^(G — G'). Detailed discussion of this numerical 
procedure is given in sec. |2.3| 

The plane-wave expansion method has been widely used in different applications, but it is not 
the purpose of this thesis to go into further detail; so far we outlined the basic numerical procedure 
of simulating a 2D PHC structure using the method. As a concluding illustration, we show in 
figure [^(b) the band diagram of a PHC of a square lattice of circular holes in a dielectric medium, 
computed with the PWE method. Fig. (a) illustrates the reciprocal lattice of such a structure, 
with the reciprocal basis vectors and the high symmetry points given for reference. 


2.2 Guided-mode expansion 

The guided mode expansion (GMF), first applied to PHCs by Andreani and Gerace [H], is, 
similarly to the plane-wave expansion method, in essence an expansion of the modes of the crystal 
over a given basis set. So far in the PWF, we only considered 2D waveguides with full translational 
symmetry along the third axis. It is only in this case that the two independent polarizations, TE 
and TM, can always be defined. In reality, however, 2D PHC-s are usually fabricated in one or 
more layers of different materials - one of the most straightforward ways to do that is by etching 
holes in an e.g. silicon slab suspended in air. Without any holes present, such a slab can be 
viewed as a waveguide due to the difference in refractive indices at the interfaces, i.e. such a slab 
supports guided modes. Thus, it seems convenient to use as an expansion basis the guided modes 
of an effective homogeneous slab of the same thickness as the PHG, but with no holes and with 
an averaged, constant dielectric permittivity 



£2(p)dp. 


( 10 ) 
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Figure 2: (a): reciprocal lattice, basis vectors, and high symmetry points of a square lattice of 
lattice spacing a, image from [3]. (b): photonic band structure of such a lattice with circular holes 
of radius 0.38a in a slab of infinite thickness and dielectric constant £2 = 9, computed with PWE. 


where A is the elementary cell area, determined based on the PHC structure in consideration. 
Below, we consider the slab to be sandwiched between two semi-infinite layers of dielectric constant 
£1 < £2 (in all the computations presented this is just air). The guided modes in such a slab have 
dispersion given by an implicit trigonometric equation, and the details for the computation of 
their energies as well as the electric and magnetic field profiles can be found in [T^] . Those modes 
are labeled by their in-plane momentum g, their mode number, a, and their polarization, TE or 
TM, i.e. we will denote the magnetic field profile of a specific guided mode as Hg ^(r) and specify 
the polarization separately (for brevity we also label fj, = (g, a)). We can then write a similar 
expansion as in eq. Q for the Bloch modes of the PHC, but expand over guided modes with 
momenta G -I- k instead of plane waves: 


Hkn(^^) — 'y ( Cr 

G,a 


(k + G, a)Hk+G,a(r). 


( 11 ) 


This summation is generally a linear superposition of both TE and TM modes, and the resulting 
Bloch mode has no well-defined polarization, i.e. all of the field components can generally be non¬ 
vanishing. However, in a system as the one we will consider, namely consisting of the same 
upper and lower claddings, those can have a quasi-TE or a quasi-TM nature, especially for low 
energies, which means that the TE (or correspondingly TM) guided modes give the predominant 


contribution to the summation in eq. (111. Moreover, the quasi-TE modes are even w.r.t. reflection 


by the z = 0 plane and are the lowest-energy ones, while the quasi-TM modes are odd w.r.t. the 
same operation. Below we present the general numerical procedure, but in all computations only 
the first, a = 1 (lowest energy) modes were used. In a symmetric structure those are TE modes, 
and thus the resulting Bloch-modes are quasi-TE. 


2.2.1 Numerical procedure 

For the computations, just as we had an eigenvalue problem in the PWE method, so do we 
obtain now, once we use eq. ([^, the expansion in eq. (0 , and the normalization condition for 
guided modes, 

I drH:(r)H^(r) = 6^,. (12) 

For the complete description with all the relevant equations of the algorithm, the reader is 
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referred to [12] . but the bottom line is that one can obtain the expansion coefficients and 

thus the Bloch-modes of the PHC structure via eq. (11), by diagonalizing the Hermitian matrix 


n 


flV — 



(VxH;(r)).(VxH,(r))dr. 


(13) 


Once the magnetic field of a mode is known, the electric field can also be computed through 
the relationship 


E(r) = X H(r). (14) 

a;£(r) 

The Fourier expansion e“^(G), as given in eq. ([^, again enters the matrix for diagonalization, 
T-Lp,u, and in the GME just as in the PWE, computing this quantity with the inverse-matrix method 
gives much better results. 

Extrinsic losses of the Bloch modes can also be treated by the GME approach, in a way 
reminiscent of the Fermi golden rule of quantum mechanics. Namely, the radiative losses can be 
represented by the imaginary part of their frequency, and this can be estimated through 9(a;) « 
3(a;^)/(2a;) and 



where A labels polarization of the radiative guided modes (TE or TM), j = 1, 3 labels whether 
those modes are leaky in the upper or the lower cladding, and pj -1- G'; is their ID photonic 

density of states at fixed wave-vector, and is the overlap between a guided and a leaky slab 
mode. 


^ ^ ^ Hi,tG'.A,,(r)) dr. (16) 

The quality factor of a resonant mode can then also be found through Q = uj/{2^{uj)) 

2.2.2 Application to waveguides and cavities 

This method can be easily used for treating defects such as waveguides or cavities in a PHC. As 
discussed before, one way to form a waveguide is to “erase” a whole row of holes, i.e. to introduce 
a 1-dimensional defect in the crystal. This introduces localized modes within the photonic band 
gap of the unperturbed crystal, in which the electric held is strongly localized to the defect region. 
In a similar way one can construct an optical cavity, by just erasing a hnite number of holes - 
one of the most studied such cavities is the L3 cavity, which is constructed by not making three 
holes in a row of a triangular-lattice-of-circular-holes PHC slab (hg. [^. In both cases, one can 
use the guided-mode expansion method to compute the PHC properties, with the difference that 
one has to consider “super-cells” instead of the normal primitive cell of the crystal in computing 
the expansion coefficients. Super-cells correspond to a collection of primitive cells, large enough 
to capture the properties of the introduced defect, and, especially for determining correctly the 
properties of cavities, large enough so that boundary effects have little to no effect. 

2.2.3 Results of computations 

We illustrate the GME method with three results of computations. First, the photonic energy 
bands of a PHC with a triangular lattice of circular holes with spacing a between neighboring holes 
was computed and can be seen in figure The parameters used were as follows: radius of the 
holes r = 0.3a, slab thickness a, and dielectric constant of the slab material £2 = 12. Furthermore, 
we limited the G-vectors in the computation up to a magnitude Gmax = 3 in units of ^. 

The method was then extended to the computation of a waveguide band diagram, the waveg¬ 
uide consisting of one row of missing holes. The super-cell considered was a square one of dimension 
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Figure 3: (a): reciprocal lattice, basis vectors, and high symmetry points of a triangular lattice 
of lattice spacing a, image from [3]. (b): photonic band structure of such a lattice with circular 
holes of radius 0.3a in a slab of thickness a and dielectric constant £2 = 12, computed with GME. 


and the result was plotted in 


figure 1^ 


The parameters used in that computation were. 


a, lO^a 

r = 0.3a, d = 0.5a, and £2 = 12. In the figure, one can see the band structure of the PHC (only 
wavevectors parallel to the direction of the waveguide are considered) and notice the two guided 
bands which lie in the band-gap of the unperturbed crystal. They correspond to a spatially-even 
and a spatially-odd band (with respect to the mirror plane). The electric field of the spatially- 
even mode with momentum vector lying on the Brillouin zone edge, i.e. fc = ^, is also plotted, to 
demonstrate the localization. 




Figure 4: (a): supercell for the PHC waveguide computation (b): band structure of the waveguide; 
the two black bands are the guided ones, (c): real-space magnetic field profile of the spatially-even 
guided mode indicated by a blue cross on the band diagram; the positions of the holes are also 
plotted 


Finally, we also illustrate an L3-cavity computation done with the GME method. Eor this to 
be correct, the supercell has to be big enough so that the field is vanishing on its edges, since 
otherwise the computation will yield the collective excitations of coupled cavities. In other words, 
what we compute is, in essence, a structure with periodically repeated supercells as the one in fig. 
each containing an L3 cavity. When they are sufficiently big, however, the results can also be 


viewed as a single-cavity computation. For 


hg. 1^ 


the supercell has size 


16a, lO^a 


with the 


other PHC slab parameters being the same as for the waveguide computation. Using eq. (151, a 
quality factor of the cavity of Q = 4659 was computed. This value is quite low, but as is discussed 
in detail e.g. in [13], it can be largely improved by slightly displacing or changing the size of 







































the hole immediately on the left and the one immediately on the right of the cavity. With such 
changes, the quality factor can go well above 10^. 
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(a) 


(b) 


Figure 5: (a): supercell used in the L3-cavity computation (b): real-space profile of the magnetic 
field of the main resonant mode inside the cavity (only the part of the supercell where the field is 
non-vanishing is shown) 


2.2.4 Discussion of PWE and GME 


The plane-wave expansion method (i.e. the form of the expansion given in eq. Q) is useful and 
very fast when computing properties of a simple structure, e.g. a regular crystal with infinite 
thickness. It can also be used to examine defects such as point defects or waveguides, but to this 
end the required computational time can significantly increase to obtain satisfactory accuracy. 
Furthermore, it is not a trivial matter to use it to model a slab of finite thickness, and even less - 
to treat radiation losses. This, however, comes quite natural in the GME method, as the guided 
modes used there are already concentrated around the core of the slab waveguide and decay in 
the claddings. This method can sometimes be more computationally demanding and a little bit 
more mathematically involved, but it has many advantages. Not only is it straightforward and 
more precise when treating a slab of finite thickness, it also performs very well when defects such 
as waveguides or cavities are introduced, and allows the computation of both mode profiles (and 
from there e.g. cavity volumes) and radiative rates (and so cavity quality factors). It is thus a 
very versatile method that can be applied in many cases and provide reliable results, while still 
keeping the numerical complexity to an acceptable level, as is expected from a basis-expansion 
method. 

Despite all its advantages, the GME is not the best tool for the purposes of this thesis - 
analyzing disorder in PHC waveguides or cavities, as, first of all, it does not allow for computations 


of big supercells, since then the diagonalization matrix in 13 becomes too large, and second of all, 
computing multiple realizations of disorder becomes cumbersome as each computation has to be 
done separately and cannot take advantage of the fact that it is only a small perturbation to a 
pre-computed structure. In sec. however, we outline the Bloch-mode expansion method, which 
remedies those problems by allowing us to treat disorder as a perturbation over a known structure, 
which could be first computed with GME. 


2.3 A numerical subtlety: “Ho’s approach” 

The “inverse matrix” method for computing £:“^(G — G') was first applied both to the plane- 
wave expansion m and to the guided-mode expansion [12] not because of rigorous mathematical 
formulation of its advantages, but rather because it seems to always give better convergence. To 
our knowledge, in the case of 2D PHCs there is, to this day, no precise mathematical theory of how 
the problem should be tackled. However, papers by Li m and Shen and He m, which concentrate 
on the ID problem, shed some light on the matter, which could better our understanding also in 
the 2D case. 
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2.3.1 Laurent’s rule vs. inverse rule for computing Fourier transforms 

The problem of computing a Fourier expansion of a periodic function is that the resulting series is 
convergent only w.r.t. the L 2 integral norm over periodic functions. In other words, the Fourier 
series of a discontinuous function is not necessarily point-wise convergent around the discontinuities 
(though the set of points where the series does not converge has to be of zero measure). The 
quantity which enters eq. (j^ is discontinuous at the hole edges. The electric and magnetic 
fields E(r) and H(r) themselves can also be either continuous or discontinuous at the boundaries 
depending on the structure and the polarization, but in the most general case they have both a 
continuous and a discontinuous component. Thus, a Fourier expansion of a product between the 
two is not necessarily point-wise convergent and thus can produce numerical artifacts. The main 
result of Li’s paper [14] is how to avoid this in a ID case: faced with a function h{x) = f{x)g{x), 
where f{x) and g{x) are in general both piece-wise smooth, periodic, and bounded, one has three 
possibilities: 1, that they have no concurrent discontinuities, 2, that they have only concurrent 
discontinuities which in the product complement each other to a continuous h(x)^ and 3, that they 
also have concurrent discontinuities which do not complement. To have best convergence in the 
numerics, it turns out that one should split the problem (starting from eq. ([^ and taking into 
account the structure and polarizations in consideration) into products of type 1 and 2 only. Then, 
for type-1 products one should apply the “Laurent rule” for computing the Fourier coefficients of 
h{x), 




M 

E 

m— — M 


fn—mQr) 


(17) 


i.e. in our case one can directly compute e ^(G) using eq. 

2, one should compute this through the “inverse matrix” rule, briefly mentioned in sec. 2.1.1 


However, for products of type 


M 

e"= E 

m— — M 



(18) 


where |/] denotes the Toeplitz matrix of Fourier coefficients of /, arranged such that the 
(n, m)-th entry of the matrix is fn-m- 

In a ID-problem, it is often possible to do such a splitting of the problem, and thus to treat the 
products “properly”. In a typical 2D case, however, for example in the case of the triangular-lattice 
of circular holes structures studied here, one would need a very complicated coordinate system with 
tangential and normal vectors to the hole edges, to be able to split the fields into fully continuous 
and fully discontinuous at the edges. This, while not necessarily impossible, would still be very 
arduous mathematically. It turns out, however, that applying only the inverse rule, sometimes 
called the “Ho method” provides very good results in many systems (and in any case, much better 
results than using only the Laurent’s rule, or a “direct computation”). As an illustration of this 
point, in fig. [^we show the guided bands of a PHC waveguide computed with the two methods, 
for Gmax = 3, 5, 7 in units of As can be seen, the convergence of the Ho-method computation 
is much better, and the same was seen for virtually all computations presented subsequently in sec. 
[^- both convergence and results reliability appeared to be much better when the inverse-matrix 
method was applied. Thus, on solely empirical grounds, everywhere in this thesis only the results 
using this method are presented. 


2.3.2 Inverting a Toeplitz-block Toeplitz matrix 

In the ID case, the matrix for inversion defined in eq. ([^ is Toeplitz, i.e. the (n, TO)-th entry is 
only a function of n — m. The numerical complexity of inverting such a matrix is greatly aided 
by an algorithm proposed by Trench [T6|, bringing it from to iV^. In the 2D case, however, 
the matrix is no longer Toeplitz because the vectors G have two components. One can, however, 
construct a Toeplitz-block Toeplitz matrix if one considers G-vectors defined on a square lattice of 
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Figure 6: The guided bands of the waveguide shown in fig. [^computed using GME where either 
a direct computation, i.e. using eq. Q (red) or the Ho method (blue) was used, for Gmax = 3 
(dashed), 5 (dashed-dotted) and 7 (solid). 


side-length 2Gmax- The vectors defined this way will have different Gx values and Ny different 
Gy values (s.t. N = NxNy), and the differences G — G' can be arranged in the following way: 


/ Gil 

Gi 2 

■ \ 

G21 

G22 • 


^Gat^I 

Gat, 2 • 

■ Gn^n^I 


where Gnm are Ny x Aty-dimensional Toeplitz matrices defined as 


(19) 


5 ^y) \^x ^x ’ ^y ^y) 


Gnm — 


(/^m /^1 ^ (/^n r^m r^2\ 

\V^ai ^y) \^x 5^2/ ^y) 




(r^n r^m /^1 r^^y\ \ 

^y i ^y ) 

(r^n r^m r^^y r^^y\ 

^x 5 '^y ^y )) 


( 20 ) 


We then need to compute the matrix e(G — G') defined over the above arrangement of G — G', 
and invert it, in order to apply the Ho method. Thus, we present the algorithm for inverting a 
Hermitian Toeplitz-block Toeplitz matrix, following the general algorithm by Wax and Kailath, 
which is itself a generalization of the algorithm for inverting Toeplitz matrices first proved 
mathematically by Trench m and given a convenient computer implementation form by Zohar, 
m- For the sake of generality, we present the inversion of a Toeplitz-block Toeplitz matrix 
composed of n -|- 1 different p x p dimensional matrices, and their Hermitian conjugates, i.e. 


( Po 
p\ 

Pi 

Pn '' 


Po 

■ Pn—1 

(21) 

\pi 

J 

Pn-1 

Po / 



where each of the block matrices poi Pi ■ • ■ Pn are Toeplitz. We start with a note on the nota¬ 
tion used: we will use the standard notations. A* and , for complex conjugation and matrix 
transpose, respectively. We will denote vectors with bold lowercase letters, and px p dimensional 
matrices will be denoted by lowercase Greek letters. Arrays of such matrices, like the matrix 
Rp^n+i given above, will be denoted by capital letters. Finally, denotes the square pxp matrix 
with units along the cross-diagonal and zeros everywhere else. This matrix plays an important 
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role here because of the persymmetric structure of a Toeplitz matrix. Note that multiplying a ma¬ 
trix on the left by is equivalent to exchanging elements which are symmetric along a bisecting 
horizontal line, while multiplying on the right exchanges the elements symmetric to a bisecting 
vertical line (and thus a matrix p is persymmetric if iTppTTp = p^). 

In general, a Toeplitz matrix is fully defined by its first row and first column. If r and c are two 
arbitrary vectors of equal dimension and with an equal first element, we will write as T = |r, cj 
the Toeplitz matrix which has r as its first row and c as its first column. Because all the matrices 
Pq ... Pn are Toeplitz, this means that we can store all the information about the matrix Rp^n+i 
in two vectors of (n -|- l)p elements, s and t, such that s contains all the firs rows of the matrices 
Pq ... Pn, while t contains all of the corresponding first columns, in the following way: 


(po)lj (Pl)l3 (Pii)l3 

S = |si . . . Sp, Sp-i -1 . . . S2p, . . . Snp+1 • ■ ■ 5(n-t-l)p} 
t — (ti . . . tp, tp-|-l . . . t2p, . . . . . . ^(n-H)p^ }- 

(Po)il (Pl)il (p„)il 

This is simply an efficient way to work with the matrices, from a computer memory point of 
view. In the description of the algorithm below we will still use the pm matrices, which are readily 
obtained: 


Pm — [[(Smp+1 • ■ ■ S(m-l-l)p)! {tmp+1 ■ ■ ■ ^(m-|-l)p)]] ^ ^ — 0 . . .Tl (23) 

The algorithm to obtain Rpl^_^_i is composed of two main parts: the first is computing the 

.. T 

p X p dimensional matrix an and the array of matrices Wn = {wi,... ojn} (corresponding to W„ 
in [HI), and the second is reconstructing recursively the full inverse matrix. We begin with the 
outline of the first part. The algorithm itself is a recursive one, where at each step am+i and 
Wm+i are calculated on the basis of am and Wm, for to = 1... n — 1. The starting values for the 
iteration are 


ai= po- pIpo Vi (24) 

Wl = {Wl} = -TTpPiCVl 

Then, for 1 < to < n — 1 we compute, with the help of the intermediate matrix Pm- 


Pm — '^pPm+1 4” ^ ( '^p^m—j + l'^pPj (2b) 

i=i 

UJj = UJj - TTpUJm-j+l{am^)'^Pm, j = 1 ... TO 
iVn-t-l — {^1 ■ ■ ■ ^m, (^m ) Pm^ 

C^m-t-1 ^m Pmi^^m ) Pm 

The full inverse matrix i?p Vi fully constructed once q:„ and Wn are obtained. In the 

case of PWE and GME computations, one needs this matrix in order to compute the matrix for 
diagonalization. However, in the Bloch-mode expansion (BME) method outlined in sec. the 
complete matrix is not needed at once, and thus as long as using computer memory efficiently is 
concerned, we would like to not store it unless we have to. Hence, we give here the procedure to 
compute it row by row, where every row is computed on the basis of the previous one. The matrix 
Rpn^i will generally not be Toeplitz, however it will still be Hermitian and persymmetric; hence 
we only need to compute one quarter of its blocks, e.g. the ones lying above both the main and 
the secondary diagonal; the rest can be expressed accordingly: 
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p—1 
-^p,n+l 


We can collect all the unique t-s in row vectors Tm = {Tm,m, Tm,m+i ■ ■ ■ Tm,n+2-m}, where m 
goes to maxrn = (n + l)/2 if n is even and to maxm = (n + 2)/2 otherwise. We can compute 
those recursively, but notice that to compute Tm+i only Tm, an and Wn are needed; thus, if 
the operations with the matrix which need to be carried out can be split into separate 

operations over the T^-s, one can optimize the usage of computational memory by only keeping 
two relatively small arrays instead of the full matrix. The rows themselves are computed by 

Ti = {Trp{a-^)^TTp, TTp{a-^fWn}, (27) 

iTnn+l)q = Tm+l,m+q = (Tm)q + ) ^m+q—1 ~ T^p^^n—m+ian U}n—m—q+2i 

l<q<n + l — 2 m, 1 < m < max^ 

For the PWE and GME computations, the main complexity still lies in diagonalizing an N x N 
matrix; in the Bloch-mode expansion method, however, using the above algorithm can improve 
the efficiency greatly, and allows for much bigger structures to be treated within reasonable com¬ 
putational time. 


t\ 2 
''”1,3 

ri .2 • • 

* '^l,n 

^l,n+l 

T2,2 •• 

'2,3 

T2,n 

'^l,n 

T2,n-1 

1 
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3 Simulating disorder with Bloch-mode expansion 

The Bloch-mode expansion (BME) [9] outlined in this section is an efficient way to treat disorder in 
PHCs, especially in structures which are too big for a supercell GME computation. Basically, the 
procedure involves computing the Bloch modes of a structure by any means (including numerical 
solution of the Maxwell equations, but here we always start from the GME-computed Bloch modes 
as they are conveniently obtained), and expanding the modes of a slightly perturbed system on the 
basis of those Bloch modes. In such a way, for example disorder in a waveguide can be treated very 
efficiently, since computing the Bloch modes of a waveguide as the one in fig. |^is relatively easy 
and requires only a small supercell, but then those modes can be used for a BME-computation of 
a waveguide containing many repetitions of the unit supercell. Furthermore, when all the Bloch 
modes are included, the computation becomes formally equivalent to a GME computation of the 
full structure; however, there is the added advantage that in many cases only some of the Bloch 
modes (e.g. the ones close in energy to the mode of interest) can be considered, and that could 
greatly diminish the computational time. Lastly, considering which of the Bloch modes overlap 
mostly with the mode of the disordered structure can provide insight into the physics of the system, 
and how it can be tuned for certain properties to be controlled. 


3.1 Theoretical basis of the BME 


To do a Bloch-mode expansion, we start from the magnetic field Bloch modes Hkn(r) of the 
regular structure, which are solutions to the equation 


V X 


1 

e(r) 


V X Hk„(r) 



Hkn(r), 


(28) 


where £(r) is the dielectric constant within the regular structure. The modes are normalized 
according to 


drH*,„,(r)Hk„(r) = 4'k^. 


(29) 


In presence of disorder, the dielectric profile e'(r) no longer has the lattice symmetries, thus 
the momentum k will no longer be conserved, so we label the eigenmodes of the system by a global 
index /3. Those eigenmodes are solutions to the equation (coming from Maxwell eqs.) 


V X 


1 


-V X H^(r) 


- ^H^(r) = 0, 


_£'(r) J c2 

and can be expanded on the basis of the Bloch modes of the regular structure: 

H; 3 (r) =^t/^(k,n)Hk„(r). 


kn 


If we insert eq. pT|) into eq. (p^, we obtain 

1 


^C//3(k,n) 


kn 


V X —V X Hk„(r) - ^Hk„(r) 


£'(r) 


= 0 . 


Using eq. (28) we can rewrite this as 


kn 


^k,n 


-I- V X ? 7 (r)V X 


Hk„(r) = 0, 


(30) 


(31) 


(32) 


(33) 


where we defined ? 7 (r) = i.e. the difference of the inverse of the regular and 

irregular dielectric profiles. If we now multiply eq. (331 by (r), integrate over r, and use the 
normalization in eq. (291, we obtain the eigenvalue problem 
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(34) 


E 

k'n' 



^kk' ^nn' + ^knk'i 


t/^(k',n') = 0, 


where the matrix Vknk'n' carries the information about disorder and is defined as 


Vk„k'„' = J drH;;,„,(r)V X (ry(r)V x Hk„(r)) (35) 

The last integral is over the volume of the considered system. We can modify it using some 
vector calculus identities and considering periodic boundary conditions, such that surface integral 
terms vanish: 


14„k'n' = I dr [(ry(r)V x Hk„(r))(V x H^,„,(r)) - V • (H;;,„,(r) x (7?(r)V x Hk„(r)))] 

= [ drr;(r)(V x Hk„(r))(V x HJ(,„,(r)) - [[ dsr,(r)HL(r) x (V x Hk„(r)) 

J J J dv 

= J drr;(r)(V x Hk„(r))(V x Hl(,„,(r)). (36) 


When ? 7 (r) is non-zero only within the slab (which is always the case if the claddings are just 
some uniform material), the above integral only needs to be evaluated there. As mentioned earlier, 
we will use GME to compute the Bloch modes, expanding the latter as in eq. (11). Thus the 
matrix elements, written using the guided modes H^(r), read 


Wn' = ^c„(m)c:,(m') /dr77(r)(V x H^(r))(V x H;,(r)) (37) 

= ^c„(/r)c;,(/r') j dr77(r)e2‘l)l^E^(r)E*,(r), 

where in the last equality we used the relation E(r) = x H(r) and the fact that for the 

guided modes e(r) is simply the average slab permittivity £ 2 - In this case summation over /i stands 
for summation over all the included reciprocal lattice vectors G and guided modes a, for fixed k 
and k'. With the expression for the electric field of a guided mode given in appendix A of [T2], and 
with the definitions of the coefficients ^ 2 ^ and l 2 ± given therein, we finally write (for a structure 
with complete symmetry w.r.t. reflection through the z = 0 plane, s.t. A 2 ^ = = B 2 ^) 


Gknk'n' = X! Cn(M)c);,(Ai')»?(g “ (38) 

2 2 

+ h-)- 


The Fourier transform 77(g — g') = £'“^(g — g') — e“^(g — gO yields much better results when 
each of the two components is computed with the Ho method, and it is here that the algorithm 
for inverting a Toeplitz-block Toeplitz matrix defined in sec. 2.3.2 comes in very handy. For large 
structures, like the ones we would like to treat with BMF, the matrices e'(g — g') and £(g — g') 
which need be inverted become very big, causing problems both from computational time and 
from computational memory perspective. The algorithm, however, solves both problems, first by 
reducing the complexity and second by removing the need to store the whole matrices at once. 

Once the eigenvalue problem in eq. (34) is solved so that the expansion coefficients 17/3(k',n') 
are known, the radiation losses of the disordered modes can also be computed in a perturbative 
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approach analogous to the Fermi golden rule in Quantum Mechanics, namely by computing the 
overlap integrals between the modes H^(r) and the leaky modes of the effective slab, 


= E E f dr^ (V X H;,(r)) • (V x 

k'n' G'ct ^ ' 

where the leaky modes (r) are given in appendix D of [T^] and are labeled by their total 
in-plane momentum g and polarization A € {TE, TM}, and j € {1, 3} distinguishes between modes 
leaking in the two cladding layers. The amplitudes 

K'", = J (V X H;,(r)) • (V x H^f/r)) (40) 

are similar to the ones computed in Appendix E of |12] , but in this case the disorder is manifest 
in the integration measure which would in general break the translational invariance of 

the system within the slab - thus, the integrals are proportional to dkk' only in the claddings. 
Therefore, for example for a structure with identical upper and lower claddings, and for TE 
polarization of both the guided and the radiative modes, and with the notation of ref. [ 12 ) . the 
amplitude is given by two terms: the claddings contribution 


.T/rad,(l) 





{Wirli+ + + VF 3 ,/i_ + A 3 ,/i+) , (41) 


and the slab contribution, 


n 


rad,(2) 
Ai'.g 


— -Cg • Eg'£ 2^2 ^(§ ~ [(fT2r + X2r) I 2 - + (^2r + lT2r) ^2+] , (42) 


where now e)"^(G — G') is the Fourier transform of the inverse of the dielectric profile within 

the cladding and e 2 ~^(g ~ s') ~ within the slab. For a more general structure and polarization 

pattern, the amplitudes can be obtained analogously, following ref. [12]. The loss rate 'yp of a 

mode is then characterized by an imaginary part of the eigenenergy ujp. This can be found by 

^2 

computing the imaginary part of given by 


rd = -3 



k.G X,j 



(43) 


and approximating 7 ^ « c^F p / {2ujp). Just as in eq. (151, pj ^g; is the ID photonic density 
of states of radiative modes, at fixed wave-vector g. 


3.2 Application to gentle confinement cavities 

As a first test of the BME method, we apply it to computing the resonant modes of three cavities 
called “Al”, “A2” and “A3” by Kuramochi et. al. [121. The cavities are based on a waveguide 
structure, in which around a certain position, the regular-structure holes are slightly displaced. 
The resulting disorder serves to localize some of the modes of the waveguide around the disorder 
position, creating cavities with very high theoretical quality factors. The cavities are illustrated 
in fig. 1 ^ (a), where besides the standard waveguide parameters we have considered so far, there 
are some additional ones. The width of the line defect in units of x/3a is labeled as tFO.98 for 
the Al cavity and tFO.9 for the A2 and A3 cavities (so essentially, the guide in fig. |^is a W1 
waveguide), while the displacement of the holes is labeled by three parameters, ds and dc 
which are however related to a single parameter x s.t. dA = x, ds = 2x/3, dc = a:/3. 
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The way we treated those structures with the BME is as follows: first, we computed the band 
structure of a regular waveguide using the GME method and a supercell of dimension 


a, 15.96^0 


for the A1 case and 


,15.8 


Vs, 


for the A2 and A3 ones. Then, using the so-computed Bloch 
modes, the procedure outlined in sec. |3.1| was applied to a waveguide composed of 16 of the 
above-mentioned unit cells, with the corresponding cavity disorder introduced in their middle. The 
parameters chosen were as in m, namely: for Al, x = 0.0214a (corresponding to the maximum 
Q found in [H]), r = 0.257a, d = 0.486a, and £2 = 11.972, and for A2 and A3, x = 0.0278a, 
r = 0.25a, d = 0.472a, and £2 = 11.972. For such parameters, the band structure is quite similar 
to the one for the guide shown in figj^ with the exception that the spatially-even guided band is 
not fully lying in the band gap (fig. ^ (a)). 

Apart from that, we restrict most of the computations to Gmax = 2( —), so that a GME 


computation with the full supercell of size 


16a, 15.96^0 


is easily carried out, for the purposes 


on comparison. In fact, when all the Bloch modes of the regular structure are included in the 
BME, the approach is formally equivalent to using GME; however, the application of the inverse 
rule separately to £'(g — g') and e(g — g') to compute r]{g — g') in eq. (39) breaks the formal 
equivalence in the finite-dimensional matrix situation (though, results are still expected to be close 
to each other). 

A natural thing to do to make a BME computation efficient is to consider only Bloch modes 
with energies close to the energy of the mode in investigation. In other words, as the waveguide 
has two guided bands in a band-gap, and the cavity mode itself will be in the band-gap (the 
main resonant mode is the spatially-even one with lowest energy), we should include the Bloch 
modes of those two guided bands, and perhaps also add the modes of the top valence and the 
bottom conduction bands. In the case of the cavities considered here, the disorder is symmetric 
w.r.t. the plane orthogonal to the slab passing through the center of the guide. Thus, the modes 
will still have the spatially-even and spatially-odd symmetry, and if we are interested in the main 
resonance mode, which is spatially-even, we do not need to include the spatially-odd guided band 
(or any other spatially-odd bands) as it will be fully decoupled. In the most basic computation, 
then, three bands were included - the spatially-even guided band, the top conduction band, and 
the lowest spatially-even valence band (which turns out to not be the lowest valence band). The 
quality factors of the cavities obtained with this computation, estimated through the loss rates 
7 / 3 , were not very satisfactory when compared to a full-supercell GME computation or the FDTD 
computation done in dg. The respective values can be seen in columns 2, 4 and 5 of table 




Figure 7: (a): band structure of the waveguide on which the Al cavity is based, the two guided 
bands are colored in black, (b): BME-computed energies of the main resonant modes of the three 
cavities normalized to the GME values for different number of bands included in the computation. 


The fact that the GME values correspond closely to the FDTD results shows that indeed the 
two methods, based on a very different foundation, agree nicely in their predictions. The fact that 
the BME computation was so far away from the GME one could only imply that more bands need 
be taken into account. Indeed, when all the 285 bands that are computed for Gmax = 2(^) were 
taken into account, the results (shown in column 3 of table came close to the GME ones, as 
expected from the mathematical equivalence of the diagonalization problems in eq. (13) and eq. 
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(39) in this case. 



3-bands BME 

285-bands BME 

supercell GME 

FDTD from [ig 

Al 

1.5 X 10^ 

7.2 X 10^ 

7.2 X 10^ 

7 X 10^ 

A2 

1.0 X 10® 

4.5 X 10® 

4.8 X 10® 

5 X 10® 

A3 

0.65 X 10® 

1.5 X 10® 

1.7 X 10® 

1 X 10® 


Table 1: Quality factors of the Al, A2 and A3 cavities with different computations 

Thus, it turns out that the high quality factor of the resonant mode is not mainly due to the 
Bloch modes lying closest in energy, but is instead due to a complex interference of all the Bloch 
modes of the system, even the ones lying very far away in the band diagram. To substantiate this 
claim, in fig. (b) we present the quality factors of the three cavities for 3, 100, 200, and 285 
bands included in the computation. It can be seen that their increase takes place gradually, and 
is not due to any single band, but is rather a complex property of the full system. However, as 
can be seen from fig. (b), the energies of the resonant modes are much better converged w.r.t. 
number of bands than the quality factors, which means that only the loss rates require many 
number of bands for an exact result. As long as energies are concerned, a 3-band computation is 
already adequate: the 3-band computed energies are correspondingly 99.94%, 99.85% and 99.90% 
of the 285-band computed ones. 

Another parameter determining the convergence of the BME and the GME alike is the number 
of G-vectors included, determined by Gmax- Therefore, in fig. [^(c) we plot the dependence of the 
3-band BME quality factor on Gmax (notice that we can relatively easily use BME to go to values 
of Gmax for which it would be almost impossible and certainly not practical to use the GME). As 
can be seen, the dependence on Gmax is much weaker than the dependence on number of bands, 
and even the Gmax = 2(— result is already realiable, as illustrated also by the fact that in the 
all-bands case it agrees with the FDTD and GME computations. 
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Figure 8: (a): the Al, A2 and A3 cavities as defined in [19]. (b): convergence w.r.t. number of 
bands of the BME-computed quality factor, normalized to the supercell GME-computed one. (c): 
dependence on Gmax of the BME-computed quality factor of the Al cavity. 


The magnetic field profiles of the main resonant mode in the three cavities are illustrated in 
fig. g for 285 bands and Gmax = 2(^). The qualitative appearance of the profiles, however, 
changes very little with the number of bands and with Gmax, in the same way that, as we already 
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mentioned, the energies change very little. One reason why so many bands are needed to get the 
correct quality factors in the case of such strong localization as is present in the considered cavities 
is that the Bloch modes are essentially extended functions. Thus, for application to cavities it 
might be wiser to expand over some more localized basis, such as the photonic Wannier functions 
[20]. However, the theoretical Q-s cited in table are never achieved in practice, due to the 
presence of random disorder. In fact, the experimentally measured quality factors cited in m 
are 800 000 for Al, 390 000 for A2, and 240 000 for A3. Conveniently enough, in the presence of 
disorder the BME seems to work well even with 3 bands, and is furthermore a convenient tool to 
quantify the effects of said disorder. 
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(a) (b) (c) 

Figure 9: Magnetic field profiles of the main resonant mode in (a) Al cavity, (b) A2 cavity, and 
(c) A3 cavity. The displaced holes are colored in the fashion of fig. @ (a). 


To illustrate this last point, disorder was included in the Al cavity computation, in the form 
of random radius fluctuations of all the wholes with mean ro = 0.257a and standard deviation 
tj = 0.0014a. The details of computing r]{g — g') in this case can be inferred from sec. 3.3.1 To 
capture better the effects of the small hole fluctuations, Gmax = 3(^) was used, and the waveguide 
length was extended to 32 elementary cells. The delicate interference which produces the very high 
theoretical Q-value of the cavity is broken by disorder, and the quality factor is brought down to 
1.07 X 10® in a computation with the 2 guided bands only (the second guided band was included 
since disorder breaks the symmetry and allows mixing). This is now quite weakly affected by 
adding more bands, with the value increasing very slightly to 1.09 x 10® in a computation with 
196 bands (all of which are spatially even in the absence of disorder, i.e. expected to couple more to 
the resonance mode than the spatially-odd ones). Thus, it is essentially only necessary to choose a 
proper disorder model to obtain numerically a certain quality factor, including the experimentally 
measured value of 800 000. And since in the presence of disorder including very few bands is 
sufficient (a claim which is also examined in the subsequent section), different disorder models can 
be handled very efficiently with a BME computation. 


3.3 Application to hole-edge disorder in waveguides 

In this section, we apply the BME to the problem of quantifying the effects of hole disorder in 
a W1 waveguide. While there is some existing literature on the matter, e.g. [H], [22] and [23], 
the approach presented here has several advantages. First, we model disorder as a general hole 
roughness given by an arbitrary deviation from a constant radius, R{4>) = Rq + SR{4>), which, to 
our knowledge, has so far never been done, while an SEM analysis of PHC structures given in 
[24] suggests that fine features of the holes are non-negligible. Second, with the help of BME we 
are able to compute very big structures, e.g. going to a waveguide of length 256 elementary cells, 
and to that we are able to average over many realizations in order to compute density of states. 
Finally, the method also allows us to compute not only the energies and radiative rates, but also 
the electromagnetic modes of the disordered guides, which had not been done prior to [3]. 
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3.3.1 Fourier transform of irregular holes 


As mentioned before, to compute the matrix elements Fkra.k'n' given in eq. (391, and specifically 
the quantity r]{g — g') = e'“^(g — g') — e“^(g — g') entering there, we need to compute e'(g — g') 
and e(g — g') as matrices and then use matrix inversion. The second of those is easily computed 
for a regular waveguide with perfectly circular holes, in the fashion of eq. ([^. Namely, let ( be 
an index running over all holes included in the supercell, with denoting their positions and 
- their radii (i.e. we allow for different radii of the holes, while still assuming they are perfectly 
spherical). Then we have 


£(g) = S25g.o + E ^ (44) 

C ^ 

where again Ji denotes Bessel function of the first kind and g is the magnitude of g. 

The quantity e^(g), on the other hand, can be very tricky to compute, and essentially depends 
on the disorder model we wish to simulate. Here, we will first assume that the disordered dielectric 
profile is constant in the z-direction, i.e. e'(r) = e{p,0), where p labels the position in the x-y 
plane. Then, we will consider a slab with holes centered at the same positions p^ (C = 1, 2,...) 
as the holes of the regular guide, but of arbitrary shape given by a contour (in polar 

coordinates), i.e. we would like to treat disorder due to irregular hole-edge. In such a setup we 
have 


-'(p) — ^2 + ~ £2)5'^(P — Pc)> 

C 


where Sq{p) is a function defined as 


0 , p>RcW 


(45) 


(46) 


To compute the Fourier transform e^(g), we start from the Fourier transform of a single function 
S'(p) as defined above, for some hole profile R{(j))'. 


r p2Tr fRW 

5(G) = / e-^^Pd^p= / d(/) 

J A Jo Jo 


(47) 


where G = (G, 6) in polar coordinates. Assume that the hole profile can be represented as 
deviating slightly from a circular shape, R{(j)) = Ro + SR{(j)). Then, 


n27T 

pRo pRQ-\-6R{(f)) 

/ 

/ +/ 

Jo 

o 

o 


= 5o(G) -f (55(G), 


(48) 


where we split the result into a circular hole expression, 5o(G) - which is well-known, and a 
correction term, (55(G). For the latter we have. 


p2Tr pRo+SR((j>) p2TT pRo+SR{<p) 

(55(G) =/ d(l) = d(f) pdpe*^'’"“(®-‘^+5), (49) 

Jo J Rq Jo J Rq 

and we make use of the following identity involving Bessel functions of the first kind: 

OO 

e*‘^“n‘^= E •/m(w)e™‘^, (50) 

m——oo 

to get 
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p27r pRQ+SR{4>) ^ 

SS{G)= / d<t^ pdp Jm{Gp)e^^(^-^^^) 

Jo Jr< 


/O JRq 

r27T oo 


m— — oo 


/o 


d</. ^ / pdpJ^iGp). 

m.^-on J Ro 


For SR small, we can expand the dp integrals to first order, 
p Ro~\~S R{4^} 


' Rq 


pdpJ^iGp) « SR{(l))RoJm{GRo). 


Since SR{4i) is a 27r - periodic function, we can take its Fourier series expansion, 


Gn 


^in4> 


n— — oo 


with reality condition G-n = G*. Plugging (52) and (53) in (51) yields 


( 51 ) 


(52) 


(53) 


^ pZTT 

dS{G) = Ro Jm{GRo) / d(/)e-™-^C'„e™'^ 

m,n— — oo ^ 

oo 

= 27 tRo Y. i^e^'"'’C^Jm{GRo) 


m— — oo 


(54) 


The Fourier expansion of the dielectric profile in presence of disorder can be directly expressed 
using this result, since 


e'(g) = £24,o + ^^5]e*®^^5c(g), (55) 

C 

where A is the area of the supercell, C runs over all the holes centered at p,^ present in the 
supercell, and 5'(;(g) = Soxis) + ^>S'c(g) is the Fourier transform corresponding to a i?f(^) = 
RqX + SR(;{(j)) profile. In the computations done here, we further simplify this result by taking 
RqX = Rq ~ constant for all holes (both in the regular and in the disordered case), i.e. only the 
hole-edge disorder SR{(j)) is different for every hole in the supercell. If we characterize this by a 
set of Fourier-expansion parameters GmX for hole number C, the expression is finally 


e'(g) = e2<5g.o + - £ 2 ) V e 


C 


,*gPC 


JiigRo) 




(56) 


3.3.2 Disorder models 

To simulate the disorder, for each hole within the supercell, a set of coefficients {Gmx} 
pseudo-randomly generated, with several requirements determining the underlying distribution. 
First, we label the standard deviation in the radius fluctuation as cr, i.e. {SR'^{(/))) = cr^, which, 
together with {G^^Cm^) ^ imposes (X)m |C'mP) = cr^- Next, an obvious choice is that the 

fluctuations increase or decrease the radius of the hole with equal probability, i.e. {5R) = 0. This 
results in (Cq) = 0. For m ^ 0, {Gm) = 0 also holds, since Cm is a complex number with random 
phase. Finally, we characterize how finely detailed the disorder is by a correlation angle 5. If we 

define exponential correlations, such that {SR{(j>)6R{(f>')) oc e“ ^ the Cm coefficients have to 
follow 
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( 57 ) 


E (|C„.P)e*”^^ = a2e-T. 

m— — oo 

In other words, the averages (ICmP) are distributed like the Fourier series coefficients of 
resulting in 

(|C„P) = ct 2 r e-Te-^-^dc/., (58) 

J — TT 

which tells us that the coefficients Cm have to be generated according to 

{Cm) = 0, {\Cm?) = , (1 - e-?(-I)"*) . (59) 

TT 1 + 

Alternatively, we could require Gaussian correlations, such that {6R{(j))6R{(j)')) oc e ^ , 

and follow the same procedure to find a different requirement for the average of the squares of the 
coefficients: 




e 25^ e 


(60) 


The result of this computation involves Dawson integrals, defined as F{z) = e e*^dt, and 


is: 


(lCmn = a^(-ir-e-^QiF 


m 


V2 V2S 


- F 


m . TT 

—=d — l—;=r- 

y/2 y/25 


(61) 


For very small values of 5, the Dawson function becomes hard to handle numerically, but the 


integrand in (60) becomes extremely well localized around 0, i.e. if we compute the integral with 


boundaries —oo to oo instead of —tt to tt, we would add only an exponentially small correction. 
The former boundary conditions are much easier to compute, and thus we can approximate 




(5—small 


(62) 


This approximation is already very good for values of (5 « but was used only for values 
i5 < 1^, where it is virtually an exact result. 

To show the difference between the two disorder models - exponential vs. Gaussian correlations 
- in fig. we plot one realization of the hole edge when either of the two correlations was used, 
in polar coordinates in (a), and the radius R as a function of the angle (j) in (b). The black line 
in fig. (b) shows the correlation length corresponding to the correlation angle S = 0.06(27r) 
which was used. The same random seed was used in both cases, which means that the same 
starting random numbers Cm with zero mean were normalized either as in eq. (59) for exponential 
correlations, or as in eq. (61) for Gaussian ones. As can be seen, the fluctuations of the radius 


with Gaussian correlations correspond much more closely to the correlation length given, while 
the exponential correlations contain fluctuations on every length scale. Furthermore, intuitively, 
Gaussian correlations seem more likely to arise in the lab because of the Central Limit Theorem. 
Thus, unless explicitly stated otherwise, everywhere below Gaussian correlations were used.. 


3.3.3 Numerical results 

As a starting point for the computations, we choose standard deviation of the radius fluctuations of 
a = 0.006a, which is consistent with (if not a bit higher than) current state-of-the art structures 
made in the lab. The correlation angle is not so easily measured and so harder to define, so 
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Figure 10: Hole profile in polar coordinates (left) and as a plot vs. (p (right) for the same ran¬ 
dom seed of Cm coefficients, but when either Gaussian-correlation scaling (green) or exponential- 
correlation scaling (blue) was used. The black line shows the corresponding correlation length. 


values of (5 = 0.06(27r) (corresponding to hole-roughness as plotted in fig. 101 and S = 0.005(27r) 
(corresponding to much finer roughness) are chosen. For completeness, computations with the 
same values, but a = 0.002a are also performed, so that the effect of the disorder “magnitude” 
a and the correlation angle S can be studied separately. To visualize the difference in correlation 
angle, in fig. El we show a realization of a hole profile with the two different values for 5, for 
a = 0.006a. 



seed of Cm coefficients, with <5 = 0.06(27r) (black) and 5 = 0.005(27r) (red). The blue dotted line 
indicates the regular hole. 


In this section, we present simulations of a W1 waveguide with 256 elementary supercells, with 
radius of the (unperturbed) holes r = 0.3a, slab thickness 0.5a, and dielectric constant of the slab 
material £2 = 12, i.e. the regular guide is the same as the one shown in fig. For the disorder 
BME computation, there are three parameters which control the precision of the computation: 
the number of bands included, the number of G-vectors defined by Cmax, and the highest number 
TUmax of moments included in the summation of eq. (|56|). The effect of all three of those is 


discussed in more detail in sec. 3.4 for a smaller waveguide (of 32 elementary supercells), with the 
assumption that a similar convergence behavior will apply. For the moment, we set Gmax = 3(^) 
and include only the two guided bands of the guide. 

To decide on a suitable rumax, we take a look at fig. |12[ where the energies and the radiative 
rates of the 20 lowest modes are plotted w.r.t. mmax, for a = 0.006a and the two different 
correlation angles S. As can be seen, for mmax = 10 the computation is fully converged in both 
cases, and even for mmax = 5 the results are converged to a satisfactory precision. Indeed, the 
radiative rates of the 6 = 0.005(27r) structure converge the slowest, but even for them 95.9% of 
all the 512 of them are within 10% of the corresponding mmax = 10 value. This precision is quite 
satisfactory both with respect to the approximations of the algorithm and with respect to the fact 
that a slightly different disorder realization can have a much bigger effect on the rates. Thus, we 
choose to stick to mmax = 5 for the other computations shown in this section. 
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Figure 12: Convergence with rrimax of the energies (left) and the radiative rates (right) of the 
lowest 20 states in the band-gap of a disordered PHC waveguide of 256 supercells, a = 0.006a 
was used everywhere, while the correlation angle is <5 = 0.06(27r) in (a) and <5 = 0.005(27r) in (b). 


This is actually a very interesting result. As can be seen in the figure, the difference between 
keeping radius fluctuations only, i.e. rrimax = 0, and including higher moments as well, can be 
quite significant: for example, all the plotted loss rates have rrimax = 0 values which are between 
10 and 50 % of the corresponding rrimax = 10 values. The overall fluctuations, however, are mostly 
due to the first 5 moments, after which the result appears to converge nicely. This shows that 
the roughest disorder in the hole edge plays the biggest part in determining the modes of the 
structure; as the features included become finer and finer, the results are affected less and less. 
This is not unexpected: due to its wave-like properties, the electromagnetic field should generally 
not be affected much by features smaller than some characteristic scale, e.g. the wavelength. 

As an illustration of the power of the BME algorithm, in fig. [^we show the magnetic field 
profiles of three different modes of a tr = 0.006a, 5 = 0.06(27r) structure. As can be seen, the BME 
method shows that the presence of disorder can induce Anderson localization of light into the 
PHC. The top image is the ground mode of the system - it has energy = 0.27177 lying below 
the band edge, and is fully localized within a region of around 10 elementary cells. The middle 
image is a mode with ^ = 0.27189, still below the band edge, but close to it. As can be seen, 
it is still localized, but within a larger area and contains two lobes - such a character is common 
for most of the states around the band edge: they exhibit two or more main lobes with decaying 
tails. Finally, for a mode with energy well within the guided band, e.g. the mode of frequency 
^ = 0.27270 plotted in the bottom of the figure, the state is extended over the full waveguide, 
even though fluctuations in intensity are of course visible because of the random disorder. 

Radiation loss rates of the modes were also computed; they are shown for all four combinations 
of disorder parameters in fig. The behavior of the loss rates for a single realization can be 
described by considering several energy regions: for lowest energies, due to the strong localization 
due to random disorder, the rates are scattered around a wide range. For slightly higher energies, 
in the region between the band edge and the light cone, the modes are no longer localized and 
show little fluctuation in 7^, up until the point when the states enter the light cone, when the 
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Figure 13: Magnetic field profile over selected portions of a W1 waveguide with 256 elementary 
cells and disorder parameters a = 0.006a and 6 = 0.06(27r). The modes correspond respectively to 
frequencies = 0.27177 (lowest frequency in the band gap), ^ = 0.27189 and ^ = 0.27270. 


radiative rates increase greatly. This is when the loss rates become of mostly intrinsic nature, 
which is also testified by the fact that they no longer depend strongly on the different disorder 
parameters. Perhaps the most striking result visible in the figure is the large effect which the 
correlation length can have on the loss rates. As can be seen, the difference in the low-modes 7 / 3 -s 
between a, S — 0.06(27r) and a <5 = 0.005(27r) structure is approximately an order of magnitude, 
and is equivalent to the difference between a tr = 0.006a and a — 0 . 002 a structure. 



o = 0.002a, 8 = 0.005 {2k) 
a = 0.006a, 5 = 0.06 {2k) 
o = 0.002a, 8 = 0.06 (271) 
o = 0.006a, 8 = 0.005 (27c) 


0 




Figure 14: Radiation loss rates of all the states inside the band gap for four different sets of 
disorder parameters, with Gaussian correlations. On the left, the bands of the regular structure 
and the light cone (black line) are given for reference. 


For the sake of comparison between exponential and Gaussian correlations, in fig. we 
present a plot of the radiative rates when the exponentially correlated disorder model of eq. (591 
was used, for the same four sets of the parameters a and <5. Two small differences between 
the results with the two models can be observed: first, the loss rates in the Gaussian case are 
slightly higher, which can be expected as for the same correlation angle, in the exponential case 
there are features much finer than the correlation length which are still included. For the same 
reason, the effect of changing <5 from 0.06(27r) to 0.005(27r) is slightly smaller in the exponential 
case, which can be inferred from the fact that the a = 0.006a, <5 = 0.005(27r) values lie slightly 
closer to the a = 0.006o,5 = 0.06(27r), as is the case for the a = 0.002a, S = 0.005(27r) and 
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cr = 0.002a, (5 = 0.06(27r) pair. 



o = 0.002a, 8 = 0.005 (27 c) 
a = 0.006a, 5 = 0.06 (2 jc) 
o = 0.002a, 8 = 0.06 {2k) 
a = 0.006a, 8 = 0.005 {2k) 
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Figure 15: Radiation loss rates of all the states inside the band gap for four different sets of 
disorder parameters, with exponential correlations. On the left, the bands of the regular structure 
and the light cone (black line) are given for reference. 


Finally, the density of states (DOS) of a structure could also be computed after averaging over 
300 disorder realizations. In panel (a) of fig. we show the DOS of the regular structure (black) 
- with a Van Hove singularity at the band edge - and of a disordered structure with a — 0.006a 
and S = 0.06(27r) (green) and <5 = 0.005(27r) (blue). For comparison, the red line shows the DOS of 
a structure in which only constant radius fluctuations were considered, with a standard deviation 
again cr = 0.006a (in other words, a case in which (Cq) = cr^, Cm/o = 0). All this is the same in 
panel (b), but for a = 0.002a. In all the cases in the presence of disorder, the Van Hove singularity 
is no longer present, in exchange for a Lifshitz tail of the DOS. The interesting result now is that 
the smearing of the energies around the band edge in the presence of disorder is strongly affected 
by the correlation length, which is obvious especially in the a = 0.006a picture (for the a = 0.002a 
picture, all the effects of disorder are scaled down and thus higher precision of the computation is 
needed for the same result to be clearly visible). All in all, we see that the correlation angle 6 of 
the hole-edge disorder can have a great effect both on the distribution of energies and of radiative 
rates. Thus, generally, ^ is a parameter which deserves consideration in the production of PHCs 
in the lab - for example, if lower loss rates are desired, decreasing S is almost as imperative as 
decreasing the overall disorder magnitude cr. 

In view of previous related work, the findings presented here are consistent with the results of 
P5] . where it was found, within perturbation theory, that local held effects due to irregular hole 
shapes blueshift and broaden the band structure. Our results extend this hnding in several ways: 
hrst, the density of states shows not only the blueshifted states, but also the redshifted ones lying 
in the Lifshitz tail. Second, the method allowed us to compute the electromagnetic prohles of the 
disordered modes. And hnally, we could also compute the radiative rates and quantify the way 
the correlation angle S affects them. Next we consider the experimental results of Le Thomas et. 
al. [^, which show that fluctuations in the hole area are the main source of degradation of the 
guided modes around the band edge (i.e. the spreading of the DOS in fig. 161. The scaling of 


disorder effects with the hole area deviations is also discussed in and [2^. While the authors 
of [25] do not observe a significant effect coming specihcally from the shape of the holes after 
the area fluctuations are taken into account, their analysis does not specifically relate to different 
correlation lengths of the disorder, but rather to some general fluctuations in the shape of the 
holes. In our model, hole area fluctuations are not separated from the shape fluctuations, but we 
see indications that different correlation lengths could have an effect, even after the area differences 
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Figure 16: (a): density of states of the regular structure (black) and of three different realizations 
of disorder with a = 0.006a. (b): same for a = 0.002a. 


are subtracted. This would make sense intuitively, since the dielectric profile is “smoothened” by 
the electromagnetic field equations, and area deviations coming from small features are “ignored”. 
Below we employ some statistics to go further into this point, specifically for the loss rates (with 
the current data, it is hard to quantify the spread of the DOS in all four cases). Within our model, 
the area of the holes has a slightly higher expectation value than the area of the regular holes, 
7 ri?g, given by 


(A) = TTcr^ + 7ri?Q. 


(63) 


Looking at the radiative rates presented in fig. 14 we see that both for 5 = 0.06(27r) and 


5 = 0.005(27r), when a is decreased 3-fold, the radiative rates change by approximately an order 
of magnitude. This is consistent with the fact that the shift in the mean value of the area in 
eq. (63) is proportional to However, the same ten-fold change is observed when changing the 


correlation angle 5 for a fixed value of ct, despite the fact that eq. (63) is independent of it. Still 


this does not mean that the area fluctuations are overall the same for different correlation angles, 
since the standard deviation a a = \/ ((^^) — (^)^) depends non-trivially on <5, due to: 


(A2)=w 



+ (Co + Ro)'^ 


(64) 


We can numerically compute this standard deviation for the disorder realizations of fig. |14[ 
to find that its change with changing either ct or <5 is almost the same. Namely, we compute 

g-A(a-=0.006a,^=0.06(271-)) ^ o c ^ <TA(g-=0.002a,^=0.06(27r)) r , r X 

crA(<T=o. 006 a, 5 = 0 . 005 ( 271 -)) ~ ~ o-a(o-=o. 002 a. 5 = 0 . 005 ( 27 r)) cases 01 cuanging 0, ana 

crA(cr=0.006a,5=0.06(277)) ^ o q ^ ita(< 7=0.006a,5=0.005(271-)) r , rases of rhanp-inp- rr Thus 

CTA(a=0.002a,5=0.06(277)) ~ c7A(a=0.002a,5=0.005(277)) CaSCS 01 CUangmg CT. iUUS, 

a change in the loss rates of a full decade between the two correlation lengths seems too much, if 
it is only due to the induced area fluctuations. Unfortunately, at this point no stronger statement 
can be made, but with a different statistical model, one can set the expectation value of eq. (63) 


to 0, and separately compare the effects of a a and 6, for a conclusive statement of weather the 
correlation length has an effect on its own. 


3.4 Discussion of the method 

The results of the previous section, while convincing and logical, are not at all conclusive because 
of the several approximations in the BME computation. Here we try to analyze those separately, 
using a 32-supercell waveguide with the same irregular-holes disorder as a toy model. In all the 
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figures here, we used a = 0.006a, <5 = 0.06(27r), but the results with the three other possibilities: 
a = 0.006a, S = 0.005(27r); a = 0.002a, <5 = 0.06(27r) and a = 0.002a, S = 0.005(27r), are presented 
in the Appendix. 

We start with the convergence of the modes with Gmax- In hg- [^(a) we plot the energies and 
radiative rates for the lowest 10 modes vs. Gmax = 3,4, 5, 6, 7 and 8(^), for rrimax = 0. As can 
be seen, increasing Gmax does not appear to change much any of the properties of the system. 
This is slightly different when we consider the same plots, but for irimax = lOj shown in Hg. 

(b), where the energies change in just the same way as in the mmax = 0 case, but the loss rates 
change more. Of course, when finer features are included in the computation, one might need a 
bit higher Fourier modes (correspondingly higher G vectors) to get convergence. Notice, however, 
that the radiative rates predominantly change in parallel. This means that one would expect the 
results with low Gmax to still be, at least qualitatively, correct. The reason why the values change 
is simply that including more terms in the expansion of the potential e^(r) corresponds in a way 
to considering a different potential, i.e. one with different local and global minima. But this is 
the same as considering a slightly different disorder realization; thus, especially for the averaging 
results taken in the computation of the DOS in fig. increasing Gmax should have little effect. 




Figure 17: (a): dependence on Gmax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a = 0.006a, 6 = 0.06(27r) and rumax = 0. (b): same 
for rUmax = 10. 


It should also be noted that the constant increase of the energies with Gmax in both cases 
is predominantly due to the increase in energies of the regular structure, rather than the BME 
computation of the irregular one. In other words, as can be seen in fig. even when using the Ho 
method, the bands are not fully convergent for Gmax = 3, and shift slightly up in energy when this 
value is increased. This is actually the main contribution to the shift of the disordered energies. 

The results of sec. |3.3.3| are very much based on the conclusion drawn from fig. that 
moments higher than m = 5 have little effect on the computation. Since we saw that increasing 
Gmax could producc an overall effect, it is imperative that we test this claim for high Gmax as 
well. To that end, in fig. 18 (a) we present the convergence with uimax of energies and rates of the 
ten lowest modes for Gmax = 3(^), while in panel (b) we show the same, but for Gmax = ’^( v)- 
As can be seen everywhere, the same result as in fig. [^is obtained: the energies and loss rates 
can change significantly with the addition of the lowest moments, but the effect is only up to 
« TO = 5; beyond that, the computation is converged. Thus, the claim that fine features of the 
hole-edge disorder do not contribute to the overall disorder effects appears to hold even when very 
high G-vectors are included. 

The last thing to verify is how the computation is affected by adding more bands. It was 
hinted in sec. |3.2| that including more bands has a huge effect on the radiative rates only, and 
only when some state resulting from complex interference is investigated, which is not the case of 
random disorder. The same result can be seen in fig. |19[ where we plot the energies and loss rates 
of a guide when different number of bands were included, for Gmax = 2(^) and mmax = 10. The 
2 -band computation involves simply the two guided bands; next the highest conduction band was 
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Figure 18: (a): dependence on rrimax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a = 0.006a, S = 0.06(27r) and Gmax = 3(^). (b): 
same for Gmax = 7(^). 


added, and afterwards the lowest valence band. Finally, computations with all bands from 1 to 
50, 1 to 100 and 1 to 185 were carried out. As can be seen, the energies stay virtually constant, 
while the radiative rates increase slightly when many bands are added. Still, the 2-band loss rates 
are « 70% of the 185-bands values, which shows again that in the case of disorder, a few bands 
are sufficient to reveal many of the properties of the system, and provide a good estimate (lower 
bound) of the radiative rates. 




Figure 19: Dependence of the energies (up) and loss rates (down) on the number of bands included 
in the BME for the 10 lowest modes of a waveguide. 
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4 Conclusions and Outlook 


The BME formalism outlined here was illustrated through application to a 2D PHC of a triangular 
lattice of circular holes in a dielectric slab. An interesting future application of the method was 
already mentioned in the end of sec. |3.3| - namely, to investigate separately the effects of hole- 
area fluctuations and hole-shape fluctuations on the modes of a structure. In addition to that, the 
method can of course equally well be applied to any 2D slab-based PHC structure, and furthermore, 
with a little extra mathematics, it can be extended to 3D structures. Here we also give two open 
problems for which it can prove very useful. 

4.1 Simulating surface roughness 

In sec. |3.3| we applied the Bloch-mode expansion formalism in order to obtain interesting results 
about how variations in the 2D-shape of holes in a PHC structure affect its electromagnetic 
modes. Another interesting question, into which the BME could give some insight, is what role 
roughness at the interfaces between the slab and the claddings plays. This has currently never 
been numerically simulated, so it is an essence still an open question. We will briefly outline how 
the BME method can be expanded to include such disorder. To that end, one can start from a 
many-layer GME computation to compute the Bloch modes: while the approach was outlined for 
three layers in [T2] (i.e. slab and two cladding layers), it is straightforward to generalize to N 
layers, as in sec. 3 of [3]. In a similar fashion, the BME formulas are straightforward to extend; 
once this is done, the simplest model one can consider is to view the regular PHC as composed 
of five layers (fig. (a)): the two claddings and the slab, plus two thin layers at the interfaces. 
The Bloch modes of such a structure can be computed with GME, and in fact, as long as the thin 
layers are from the same material as the slab, a 3-layer GME calculation is sufficient. We can 
then simulate surface roughness by leaving the slab (layer 2) and claddings unperturbed, while 
changing the profile of the thin layers (1 and 3) to some arbitrary shape (fig. (b)). As long as 
this change does not depend on z, it can be conveniently treated in a BME formalism, with the 
matrix Ekrak'n' determined by the disordered profiles of the two thin layers, ei(r) and e 3 (r). The 
proposed model, while very simplified, could provide plenty of insight into the effect of surface 
disorder 

Upper 

cladding 

Layer 3 -C 

Layer 2 - 

Layer 1 -c 

Lower 

cladding 

(a) (b) 

Figure 20: (a): cross-section of the PHC without disorder: it can be viewed as composed of 5 
layers instead of 3. (b): adding random fluctuations (constant in the vertical direction) to layers 
1 and 3 simulates roughness at the interfaces. 


on the PHC modes. 
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cladding 



4.2 Quantum dots embedded in a PHC structure 

Quantum dots embedded in PHC structures promise to open up almost unlimited possibilities 
for applications starting from simple light emission and/or lasing to implementation of qubits 
and qubit gates for quantum information processing. With the ever-advancing sophistication of 
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our ability to produce such integrated structures in the lab [27j . it is imperative that numerical 
simulations follow suit. The electromagnetic modes of a PHC structure which can be computed 
in a BME formalism can be the starting point for such simulations. For example, in a Maxwell- 
Schrodinger framework as the one developed and applied in |28j . one models the presence of the 
dots by adding a (generally nonlocal) linear susceptibility tensor y(r, r', cu), which can be computed 
to a different degree of sophistication based on the quantum-dot wavefunction models one wishes 
to include. This, together with the electromagnetic mode of the background medium is sufficient to 
explore the radiative properties both of a single quantum dot (energy-shift and radiative lifetime) 
and of a collection of dots (radiative coupling and collective modes). The fact that the BME makes 
it easy to compute small variations to a starting structure, then, would also make it possible to 
optimize structures for the tuning of the quantum dot properties. 

4.3 Conclusion 

In this thesis, we summarized some basis-expansion methods for numerically simulating 2D Pho¬ 
tonic Crystals, namely plane-wave expansion and guided-mode expansion for simulating regular 
periodic structures, and Bloch-mode expansion for the addition of disorder. The main results of 
this work are based on the Bloch-mode expansion, a method which allows for very large structures 
(compared to the limitations of e.g. FDTD or supercell-GME) to be simulated, and yet can be 
used to accurately compute not only the frequencies and loss rates of the modes, but also their 
real-space profiles. In sec. we first tested the method against some gentle-conhnement cavi¬ 
ties for which results already appear in the literature, as well as against guided-mode expansion 
with a big supercell, and saw that in the case of strongly-conhned modes, the BME produces 
quantitatively reliable quality factors only when many bands are included, even though even a 
few bands can give good information about the mode prohles and their energies. In the presence 
of disorder, however, BME appears to be fully quantitatively reliable even for a low number of 
bands, as illustrated by the small change in quality factors of the cavities between a 2-band and 
a 195-band computation with disorder. 

A similar observation was made for the case of disorder in W1 waveguides, and a 2-band BME 
was applied to the problem of quantifying the effect of hole-edge disorder on the modes of the 
waveguide. This is still an open question, and the results obtained are the main innovation of this 
work. Namely, we found that if we split the hole prohle into Fourier components, only the hrst 
« 5 moments can have a strong influence on both the energies and the loss rates of a guided mode, 
and for rrimax = 10 the computations were always fully converged. The explanation to this lies 
most probably in the fact that sufficiently fine structures become “invisible” for an electromagnetic 
wave of a certain wavelength, i.e. the underlying Maxwell equations smoothen the potential so 
that only the big features matter. This resulted in the fact that for the same magnitude of 
disorder, tr = 0.006a or cr = 0.002a, a difference of correlation angle between d = 0.06(27r) and 
5 = 0.005(27r) can produce an order of magnitude of difference in the loss rates of the guided 
modes. Furthermore, for higher correlation angle, we saw a stronger smearing of the DOS of the 
waveguide around the guided-band edge, which suggests stronger localization of the modes and 
so degradation of the dispersive slow-light regime. It will be interesting to extend this work by 
rehning the statistical model of disorder and considering more sets of the disorder parameters. 
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A Results as in sec. \3A[ but for different disorder param¬ 
eters 


and rumax, 
a = 0.006a, S = 


22 and 251, and a = 0.002a, i5 = 


In this appendix we present data related to the convergence of the BME vs. Gmax 
similar to the data presented in sec. |3.4| but for the three remaining cases: 

0.005(27r) (fig. J2^and@, cr = 0.002a, ,5 = 0.06(27r) (f “ “ 

0.005(27r) (fig. |23| and 261. The discussion from sec. |3.4| applies equally well (and sometimes 
better) in all cases. Most importantly, the fact that the moments up to « m = 5 have the main 
contribution to energies and loss rates is visible in all cases shown in fig. and 1^ 



0.28 

^ 0-275 



X 10' 



Figure 21: (a): dependence on G^ax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a = 0.006a, 6 = 0.005(27r) and nimax = 0. (b): 
same for nimax = 10. 




Figure 22: (a): dependence on Gmax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a = 0.002a, 6 = 0.06(27r) and nimax = 0. (b): same 
for nimax = 10. 
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Figure 23: (a): dependence on Gmax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a = 0.002a, 5 = 0.005(27r) and rrimax = 0. (b): 
same for rrimax = 10. 



Figure 24: (a): dependence on rrimax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a — 0.006a, S — 0.005(27r) and Gmax = 3. (b): 
same for Gmax = 7. 



Figure 25: (a): dependence on rrimax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with cr = 0.002a, 6 = 0.06(27r) and Gmax = 3. (b): same 
for Gmax — 7. 
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Figure 26: (a): dependence on rrimax of the energies (up) and radiative rates (down) of the 10 
lowest modes of a 32-supercell W1 guide with a — 0.002a, i5 = 0.005(27r) and Gmax = 3. (b): 
same for G^ax = 7. 
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